Información sobre el dataset

Vamos a utilizar el dataset GSE7614. Este dataset corresponde al estudio Phenotypic and transcriptional response to selection for alcohol sensitivity in Drosophila melanogaster.

Información del dataset
Organismo Drosophila Melanogaster
Número de muestras 24 (15 individuos por muestra)
Grupos Macho Control (x2 Rep. Biológicas x2 Líneas)
Hembra Control (x2 Rep. Biológicas x2 Líneas)
Macho Sensible (x2 Rep. Biológicas x2 Líneas)
Hembra Sensible (x2 Rep. Biológicas x2 Líneas)
Macho Resistente (x2 Rep. Biológicas x2 Líneas)
Hembra Resistente (x2 Rep. Biológicas x2 Líneas)
Número de características 18952
Plataforma [Drosophila_2] Affymetrix Drosophila Genome 2.0 Array
GPL1322
Tipo de muestra in situ oligonucleotide

Cargado de datos crudos

En este caso hemos descargado todos los archivos .CEL del dataset GSE7614. También hemos consultado a que grupo pertenece cada una de las 24 muestras que contiene el dataset. Hemos preparado también un archivo que contiene toda la información referente a los grupos y las muestras y lo hemos cargado en memoria.

short_name sigue la convención c(ontrol)/r(esistant)/s(ensible) m(ale)/f(emale). Los números separan las diferentes samples y las grupos de los que salieron (ver información sobre el dataset).

Luego hemos generado un ExpressionFeatureSet juntando los archivos .CEL y annotatedDataFrame que contiene el phenoData

Tenemos en total 535824 filas/sondas y 24 muestras.

Control de calidad de los datos crudos

Para realizar un control de la calidad de los datos crudos vamos utilizar la libreria Array Quality Metrics y affyPLM.

Lo haremos de dos maneras diferentes, primero sin transformar los datos con logaritmo de la intensidad y luego haciéndolo.

Sin transformar

Salida Array Quality Metrics sin transformar

Salida Array Quality Metrics sin transformar

  • En el gráfico A podemos observar un resumen de todas las muestras y el resultado obtenido del control de outliers. Vemos que según el algoritmo las muestras 1 y 17 podrían presentar algún problema por ello las tendremos en especial consideración cuando observemos los gráficos. En este caso no disponemos de la posibilidad de hacerlo, pero podríamos comprobar las imágenes de los microarrays para ver si ha ocurrido alǵun problema durante su procesado.

  • En el gráfico B observamos una agrupación por distancias de los microarrays, vemos que hay una separación clara según el sexo de la muestra, a excepción de las muestras 1 y 17 que muestran una distancia un poco atípica se encuentran más distantes del resto de muestras de machos.

  • En el gráfico C observamos los boxplots de las intensidades de cada uno de los arrays. También aparece marcado el array 17 quizá más disperso que el resto, también aparece marcado el 10, sin embargo no parece que sean muestras particularmente aberrantes.

  • En el gráfico D podemos observar un análisis por componentes principales, veremos debajo uno como el que se ha realizado durante las prácticas de la asignatura, en este caso vemos una separación clara entre los grupos macho y hembra de las muestras lo cual es esperable y deseable.

  • En el gráfico E, vemos que las muestra siguen un mismo patrón, aunque la escala hace un poco complicada esta observación, probaremos luego con una función de affysimple.

  • En el gráfico F, el valor esperable es que la mediana sea constante, sin embargo aparece una pequeña joroba en el lado derecho, pero seguramente esté relacionado con una saturación de la intensidades como indica el propio paquete.

  • En el gráfico G, podemos observar la comparación de cada uno de los arrays respecto a un pseudo-array que tiene el valor mediano a través del resto de arrays. En este caso lo esperable es que los valores se concentren a lo largo del eje M=0. Sin embargo al no estar transformados los valores aparece agrupados en el margen izquierdo haciendo difícil la visualización.

Vamos ahora a ver un PCA más completo

Vemos una agrupación clara por la condición sexo, no tanto por las condiciones experimentales.

Log-Transform

Array Quality Metrics Output

Array Quality Metrics Output

Al estar transformados con el logaritmo la forma de las gráficas cambia ya que hemos estirado la distribución. El significado de los gráficos es el mismo que en el apartado anterior.

El array 17 sigue apareciendo como un outlier, sin embargo en ninguno de los plots aparece particularmente dramático por tanto lo conservaremos.

affyPLM

Vamos ahora a utilizar la librería affyPLM vamos a ajustar un modelo a nivel de sonda y vamos a comprobar las gráficas RLE y NUSE.

En el caso del RLE vamos a encontrar una medida estandarizada de la expresión, que ha de ser similar a través de arrays, y en el caso de NUSE observaremos los residuos normalizados del modelo cuya distribución ha de ser similar a través de arrays.

No parece que ninguno de los arrays se desvíe particularmente en ninguna de las dos gráficas (quizá la c_1_m_1), si que vemos que hay una diferencia clara entre sexos. Podríamos preguntarnos en este caso si se ha cometido algún error y los machos perteneces a un batch y las hembras a otros, es decir hemos generado algún tipo de confound.

Conclusiones de calidad

  • Conservamos todos los arrays
  • Podría ser útil revisar el 17
  • Podría ser útil consultar si ha habido alguna diferencia sistemática entre los arrays de machos y hembras, batch, fecha, operario,…

Normalización

Procedemos ahora a normalizar los datos con el método RMA.

Durante este paso el algoritmo:

  • Eliminará el ruido de fondo
  • Normalizará la señal entre arrays, nunca entre condiciones experimentales.
  • Colapsará todos los valores de cada grupo de sondas en un único valor por gen.

Una vez normalizado vamos a volver a realizar un análisis de calidad de los datos

Array quality metrics output después de la normalización

Array Quality Metrics Output

Array Quality Metrics Output

Los gráficos son los mismos que en apartados anteriores, ahora vemos una separación mucho más clara en el array de distancias. El array número 17 ahora es más similar al resto.

Solo existen dos detalles destacables, uno es la diferencia entre patrones de intensidad en la figura E entre las muestras machos y hembras, y la joroba que aparece en la figura F, cuando la mediana debería de ser plana, quizá sea debida a la diferencia entre sexos.

Conclusiones normalización

Los datos parecen correctos, por tanto no es necesaria ninguna acción.

Filtrado

Vamos a filtrar las probes que no tengan entrada en y las que corresponden a las sondas de control de Affymetrix. Vamos a hacerlo de manera manual utilizando la base de datos correspondiente al chip: Affymetrix Drosophila Genome 2.0 Array, en el paquete ‘drosophila2.db’. Lo que haremos será eliminar todas las features cuyo entrada EntrezID sea NA.

También vamos a eliminar algunas sondas que tienen varios genes asignados, o los genes que tienen varias sondas asignadas. Puede que exista alguna manera más correcta de realizar este paso, pero por sencillez las vamos a eliminar.

De esta manera pasamos de tener 18952 features a 11294. Lo cual nos ayudará, por ejemplo, reduciendo el impacto de la correción por múltiples comparaciones.

Modelo y contrastes

Estamos ante un diseño de 3x2, en este caso vamos a simular que estamos interesados en las comparaciones:

  • Main Effect Macho vs. Hembra
  • Single Effect Resistente vs Sensible
  • Interacción Sexo vs (Resistente vs Sensible)

Vamos a ajustar un modelo lineal con Limma y la siguiente matriz de diseño

Horizontal axis indicates the term of the design matrix. Vertical axis indicate the short name of the sample

Horizontal axis indicates the term of the design matrix. Vertical axis indicate the short name of the sample

Una vez ajustado el modelo calculamos los siguientes contrastes

Anotación

Anotamos los datos con el paquete drosophila2.db.

Resultados

Genes diferencialmente expresados entre condiciones

Las tablas permiten ordenar y filtrar los genes en función de diversos parámetros. Solamente se incluyen los 1000 genes más diferencialmente expresados por cada comparación.

Corrección por multiples comparaciones False Discovery Rate Benajin y Hochsberg CITA!

Macho vs. Hembra

Resistente vs. Sensible

Interacción Sexo y (Resistente vs. Sensible)

Volcano plot

Aquí presentamos un volcano plot con algunos extras sobre el original, con los p-valores ajustados y sin ajustar. Además se incluye una tabla con el porcentaje de genes diferencialmente expresados en cada uno de los contrastes. Se incluye el corte en p<.05 para valores ajustados y sin ajustar.

Volcano plot modificado. Las líneas horizontales indican p-valor=.05. Las líneas  verticales indican el un fold-change de 2.

Volcano plot modificado. Las líneas horizontales indican p-valor=.05. Las líneas verticales indican el un fold-change de 2.

Porcentaje de genes significativos
contrast Sin corregir p<.05 Corregido FDR p<.05
int_sex_resvssens 5.2% 0.13%
male_vs_female 86% 86%
res_vs_sens 19% 4.1%
a Redondeado a dos decimales

Diagrama de Venn

Los resultados que aparecen aquí están corregidos al p<.05, con el p-valor corregido a través de múltiples comparaciones con FDR (Bejamini Hochberg) y se han considerado los contrastes de manera separada para esta correción.

Número de genes diferencialmente expresados por contraste
res_vs_sens male_vs_female int_sex_resvssens
Down 198 5825 9
NotSig 10833 1619 11279
Up 263 3850 6
Diagrama de Venn que muestra el número de contrastes significativos encada uno de los contrastes y sus intersecciones.

Diagrama de Venn que muestra el número de contrastes significativos encada uno de los contrastes y sus intersecciones.

Tabla de anotaciones

Se ha generado (aunque no sea accesible desde el documento, si que está en el repositorio de git) con todos los genes anotados que han resultado signficativos en cada uno de los contrastes.

Archivo: ‘annotations.html’

TODO ANOTAR LOS GENES CON MAYOR FOLD CHANGE

Perfiles de expresión

Análisis de enriquecimiento

  • Para estos análisis hemos seleccionado todos aquellos genes cuyo p-valor corregido con FDR es menor a .1.
    • Hemos realizado un análisis por contraste
    • Hemos utilizado las bases de datos KEGG y GO
    • En el caso de KEGG hemos corregido siempre con un p<.1
    • En el caso de GO hemos corregido siempre con un p<.1 y un q<.05
    • Hemos rebajado ligeramente el umbral de significación para poder mostrar los gráficos y los tipos de análisis disponibles, de otra manera esta sección quedaría vacía.
      • En un análisis real mantendríamos el umbral de signifiación fijado a priori, o marcaríamos el análisis como exploratorios.
    • Hemos realizado 10000 permutaciones en los GSE

    Enlaces a los tipos de gráficos utilizados:

Gene Ontology

Explicación de para que vale la Gene Ontology

Over-representation analysis

Solo cálcula si es por azar que los genes sobre representados aparecen más en un pathway que en otro, con un fisher o un binomial, por azar deberían de aparecer, pero en realidad aparecen.

Macho vs. Hembra

En este caso aparecen muchos genes relacionados con cada ruta y la visualización se hace un poco complicada. Un upsetplot hubiese sido más útil, pero en este caso el paquete daba un error de compatibilidad que no he sido capaz de solucionar.

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resistente vs. Sensible

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Interacción Sexo y (Resistente vs. Sensible)

## [1] NA

Gene Set enrichment analysis

Aquí el approach es un poco distinto porque tenemos en cuenta si se están up o down regulated todos los metabolitos a cañón ahí dentro.

Macho vs. Hembra

Resultados de GSE para GO

Resultados de GSE para GO

Resultados de GSE para GO

Resultados de GSE para GO

Resistente vs. Sensible

## [1] NA

Interacción Sexo y (Resistente vs. Sensible)

## [1] NA

KEGG

Over-representation analysis

Contrast Rutas sobrerrepresentadas
Macho vs. Hembra 2
Resistente vs. Sensible 0
Interacción Sexo y (Resistente vs. Sensible) 4

Macho vs. Hembra

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resistente vs. Sensible

## [1] NA

Interacción Sexo y (Resistente vs. Sensible)

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Gene Set enrichment analysis

En los casos en que hemos encontrado rutas significativas, hemos utilizado el paquete pathway para descargar la ruta de la base de datos de KEGG e incluirla en el documento. Así los genes significativos se pueden observar dentro del contexto de la ruta.

Contrast Rutas sobrerrepresentadas
Macho vs. Hembra 0
Resistente vs. Sensible 2
Interacción Sexo y (Resistente vs. Sensible) 0

Macho vs. Hembra

## [1] NA

Resistente vs. Sensible

## Picking joint bandwidth of 0.926
Resultados de GSE para KEGG

Resultados de GSE para KEGG

Resultados de GSE para KEGG

Resultados de GSE para KEGG

Metabolism of xenobiotics by cytochrome P450 Drug metabolism - cytochrome P450

Interacción Sexo y (Resistente vs. Sensible)

## [1] NA